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Abstract 

The accuracy and reliability of gamma-ray bursts (GRBs) as distance indicators 
are strongly restricted by their systematic errors which are larger than statistical er- 
rors. These systematic errors might come from either intrinsic variations of GRBs, 
or systematic errors in observations. In this paper, we consider the possible origins 
of systematic errors in the following observables, (i) the spectral peak energies (E p ) 
estimated by Cut-off power law (CPL) function, (ii) the peak luminosities (L p ) esti- 
mated by 1 second in observer time. Removing or correcting them, we reveal the true 
intrinsic variation of the E v -T^-L v relation of GRBs. Here Ti is the third parameter 
of GRBs defined as Tl = Ei SO / L p . Not only the time resolution of L p is converted from 
observer time to GRB rest frame time, the time resolution with the largest likelihood 
is sought for. After removing obvious origin of systematic errors in observation men- 
tioned above, there seems to be still remain some outliers. For this reason, we take 
account another origin of the systematic error as below, (iii) the contamination of 
short GRBs or other populations. To estimate the best fit parameters of the E p -Tl- 
L p relations from data including outliers, we develop a new method which combine 
robust regression and an outlier identification technique. Using our new method for 
18 GRBs with <je p /E p < 0.1, we detect 6 outliers and find the E v -T^-L v relation 
become the tightest around 3 second. 

Key words: gamma rays: bursts — gamma rays: observations — gamma rays: 
cosmology 
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1. Introduction 



There are some correlations between the spectral peak energies (E p ) and the brightness 
of gamma ray bursts (GRBs), as well as correlations between temporal properties of bursts 
(variability, spectral lag, etc) and their brightness (Fenimore & Ramirez-Ruiz 2000; Reichart 
et al. 2001; Norris et al. 2000). The first correlation found is the E^-isotropic equivalent energy 
(E iso ) correlation found by Amati et al. (2002). This correlation is confirmed and extended to 
both higher and lower energy by many satellite teams (Amati 2006; Sakamoto et al. 2004; Amati 
et al. 2009; Krimm et al. 2009). Although there are some paper against the E p -E iso relation, 
all of them use poor data without redshift and/or spectropic E p (Butler et al. 2007; Ghirlanda 
et al. 2008; Nava et al. 2008). It is difficult to know whether these results come from the 
intrinsic property of GRBs or the systematic effect by using poor data. 

Another correlation is the -E^-jet collimated energy (i? 7 ) correlation found by (Ghirlanda 
et al. 2004). This is one of the tightest correlation and expected to be used as distance indicator 
toward high redshift universe (Ghirlanda et al. 2006), but the need for jet break time make 
difficult to increase their sample. More importantly, the complexity of early afterglow observed 
by Swift makes it difficult to identify the jet break time. Thus, it seems to be not so effective 
to use the E p -E^ correlation and also Ep-T^^-Eiso correlation (Liang & Zhang 2005) for 
cosmology. 

A correlation between E p and their 1-second peak luminosity (L p ) was also found by 
(Yonetoku et al. 2004). Firmani et al. (2006) proposed that taking the high signal time (T0.45) 
as third parameter improves the E p -L p correlation. Further, Rossi et al. 2008 showed that 
the Ep-To^s-Lpcorrelation becomes as tight as E p -E- lso correlation, using 41 GRBs observed by 
BeppoSAX and Swift . Recently, Tsutsui et al. (2009) found similar but different correlation 
between E p , luminosity time (T L = E iso /L p ) and L p . Because these relations use only prompt 
emission parameters, they might become useful distance indicators for high redshift universe. 
However, there are some researches against these correlations. Collazzi & Schaefer 2008 studied 
whether adding any timescale improves the E p -L p correlation (not only T 45 , but also T 90 , T L , 
and so on.), but they found that there is no timescale which improves the relation. From the 
result they insist that the E p -T 0A5 -L p (E p -Tl-L p ) relation is equivalent to the E p -L p relation. 
However there seems to be some room for further investigation. 

First of all, we must point out that there are many reasons which cause the scatter 
around the relation in addition to the intrinsic dispersion and the measurement uncertainties. 
As discussed in some previous works (Kaneko et al. 2006; Krimm et al. 2009; Shahmoradi & 
Nemiroff 2010), the peak energies estimated by the Cut-off power law (GPL) model become sys- 
tematically higher than the peak energies estimated by the Band model. Besides, the timescales 
of GRBs must be defined in GRB rest frame, because fixed observed timescales correspond to 
different rest-frame timescales for GRBs at different redshifts. Then the parameters defined by 
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observer time also would cause extra scatter and redshift dependence of these relations (Tsutsui 
et al. 2008; Yonetoku et al. 2010). For these reasons, we use only the GRBs whose peak energies 
are estimated by the Band model, and peak luminosities computed using a fixed timescale in 
GRB rest frame in this paper. 

Even if we take these factors into account, another problem still remains: how many 
populations of GRBs exist? Clearly there are some GRBs which are far from most of GRBs 
in the relations, e.g. GRB980425, and some short GRBs. Furthermore, there might be other 
unknown populations. Therefore, when we derive correlations, applying ordinary regression 
to the data including different populations would result in misleading results. Thus, we must 
separate different populations of GRBs in our regression analysis, or must use robust statistics. 
As far as we know, there are no reference to deal with data which potentially contains both 
intrinsic dispersion and/or multiple populations. In this paper, we develop a new method to 
do this. 

The structure of this paper is as follows. First we describe our database of 86 GRBs with 
known redshift, public light curve and their spectrum parameters (§-2). In §-3, we describe 
the new method we developed. Using our new database and method, we estimate the best- 
fit relation and true intrinsic dispersion in the E p -Ti J -L p plane (§-4). In §-5, we present and 
discuss our results. 

2. Data Description 

In Yonetoku et al. (2010), we constructed a database selecting 109 GRBs from GCN 
Circular Archive (Barthelmy 1997) and GRBlog (Quimby et al. 2004). In this section, we 
briefly describe our database. 

In many cases, the prompt gamma-ray spectrum is well fitted with the spectral model of 
the exponentially-connected broken power-law function suggested by Band et al. (1993). This 
Band function has four parameters, the low-energy photon index a, the high-energy photon 
index (3, the spectral break energy E and the normalization. The peak energy (E p ), at which 
the flux is maximum in the vF v spectrum, can be calculated as E p = (2 + o)Eq. 

However, for some GRBs, the photon index (mostly 0) cannot be determined due to the 
limited energy range of the detector and/or the lack of the number of photons (Pendleton et al. 
1997). When the observation of high-energy range is not enough, the spectrum is sometimes 
fitted with the cut-off power-law (CPL) function. This CPL function is very similar to the 
Band function, but the high energy end is exponentially cut-off. This function is composed 
by three parameters, the low-energy photon index a, the spectral break energy E and the 
normalization. The peak energy can be also expressed as E p = (2 + a)Eo similar to the Band 
function. Note that, for an observed GRB spectrum, even if the reduced x 2 value of the CPL 
model is smaller than that of the Band function, it is difficult to say whether this model reflects 
the intrinsic property of the GRB or it is just due to the poor statistics in the high-energy 
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range. 

In Yonetoku et al. (2010), we calculated the bolometric energy and the peak luminosity 
in the energy range 1-10,000 keV in the rest frame of each GRB by extending the observed 
spectrum. Here, it should be noted that the integration was performed assuming the Band 
function even for GRBs whose spectra were not fitted by the Band function and the high 
energy photon index were not reported. In these cases we assumed the typical values a = — 1 
and j3 = —2.25 to calculate the bolometric fluence (Sboi) and the bolometric peak flux F p boi- 
These values are suggested by BATSE observations (Preece et al. 2000) and also supported 
by Fermi observations of GRBs up to 100 GeV energy range. Zhang et al. (2010) confirmed 
that the time- resolved spectra of 14 out of 17 GRBs are best modeled with the Band function 
over the entire Fermi spectral range. Then the bolometric isotropic energy (E lso ) and the 1- 
second peak luminosity (L p ) can be simply calculated as E iso = And 2 L St, Q \/ {I + z) (erg), and 
L p = And^Fp^oi (erg s~ x ). Here, di is the luminosity distance calculated for the flat A-CDM 
universe with the cosmological parameters of (fi m ,f2A) = (0.3,0.7) and the Hubble parameter 
of Hq = 70 km s _1 Mpc _1 . Further we define the luminosity time as the third parameter of 
GRB prompt emission as Tl = E- lso / L p . The error of the luminosity time is estimated by using 
error propagation equation. We can neglect the crossterm between L p and E iso because of the 
independence of the E p -E- lso and E p -L p relation shown in (Tsutsui et al. 2009). 

Thus, for GRB whose observed photon number is small, there are two possible systematic 
effects. One comes from the fact that the peak energy E p is determined by fitting the spectrum 
with either the Band function or GPL function. As Kaneko et al. (2006) pointed out that the 
CPL function tends to overestimate E p compared to the Band function. This would induce a 
systematic error in the correlations related to E p . On the other hand, although L p and Ei SO are 
determined in a single straightforward way, the photon indices are set to the typical values if 
the number of detected photons is small. This would also cause a systematic error. 

To estimate the impact of systematic errors, we separate GRBs into three group. One 
is the platinum data set which consists of GRBs whose spectrum is well observed and fitted 
by the Band function with small E p error, cte p /E p < 0.1. The second, gold data set, is defined 
by GRBs fitted by the Band function with relatively large error, (Te p /E p > 0.1. Finally bronze 
data set consists of all other GRBs. 

In the previous works, the peak luminosity was culculated with 1-second peak flux in 
observer frame. This means that the time scale of the peak luminosity is different from event by 
event in GRB frame because of the different redshift. So the time scale should be defined in GRB 
rest frame. Then we must convert the peak luminosity reported in Yonetoku et al. (2010) to the 
r-second peak luminosity (L p (t)) in GRB rest frame. Let us explain the conversion method. 
The archived lightcurves of CGRO-BATSE and Konus-Wind are provided with 64 msec time 
resolution, and we can also create the lightcurve of Swift-BAT with the same time resolution. 
We used only 86 GRBs observed by these three satellites in this paper, although there are 109 
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GRBs in the database of Yonetoku et al. (2010). This is because there are no lightcurves with 
64 msec resolution for the remaining 23 GRBs mainly observed by BeppoSAX or HETE-II. 
When we consider the r-second peak luminosity in GRB frame, the time scale of observed 
peak flux is equivalent to r(l +z)-second because of the cosmological time dilation. Then we 
performed the re-binning with the number of N(t) = r(l + z)/0.064 bins (round off to the 
nearest whole number), and estimate the observed peak photon flux. Here, we have N(t) 
degrees of freedom to choose the start point of re-binning. Therefore we searched all patterns 
of re-binning to find the brightest peak photon flux P(t(1 + z)). We convert the bolometric 
peak flux F p>bo i(r(l + z)) with 

F pM (r(l + z)) = W1.024) x P p^y W 

Here we refer to the 1-second bolometric peak flux F p ^ i (1.024) for each event summarized in 
Yonetoku et al. (2010). The peak luminosity calculated by L p (t) =47rd'^F Pi b i(T(l + z)) becomes 
the same time interval of r-second in GRB frame. Therefore L p (r) may be more appropriate 
than the previous definition by Yonetoku et al. (2010) to discuss the E'p-brightness correlations 
of GRBs. In the following sections, we use newly estimated L p (t) as the bolometric peak 
luminosity L p . 

In the whole analysis of this paper, we do not use GRB980425 which is a famous outlier 
for both E p -L p and E p -E lso relations. The data are summarized in table. 1-3. 

3. Correlation analysis with outliers and multiple populations 

Ordinary regression assumes that the scatter of data around the relation follows a 
Gaussian distribution. However, in many cases of astrophysics, there are some obstacles which 
lead to misleading results. For example, as in the case of Cepheid variables, there might be sev- 
eral populations whose relations are very different from each other. Also there might be some 
experimental mistakes which result in outliers. In these cases, we must separate each popu- 
lation or identify and remove outliers, because they would dominate the chi square value and 
bias the result of regression. In the history of astronomy, such kind of mistakes was frequently 
seen. For example, it is very famous that Hubble combined different types of Cepheid variables, 
without knowing it, to estimate the distance to galaxies. Consequently he overestimated the 
expansion rate of the universe by a factor of two. Likewise, there are some clear outliers in the 
correlations of GRBs, such as GRB980425 and some short GRBs. Thus, we should remove the 
outliers to obtain the true correlations and estimate the intrinsic dispersion reasonably. 

If we knew a criterion to distinguish different populations and the value of the intrinsic 
dispersion of the relation in advance, there are some reliable ways to eliminate outliers (Kowalski 
et al. 2008). However, we cannot find any reference to do this without any prior knowledge. 
Because the number of GRB events with small observational errors is currently very small, it 
is important to eliminate outliers in a systematic manner rather than in an ad hoc manner as 
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is often seen. Here we develop a statistically reliable way to derive a correlation and estimate 
the intrinsic dispersion. 

The basic idea is to combine a robust regression with outlier detection (Hampel et al. 
2005). To identify outliers, we have to evaluate the residuals of samples from the best-fit 
relation. However, if we perform ordinary chi square regression, outliers can influence the result 
of regression itself so much that we cannot measure the residuals correctly. Thus, we should 
adopt a more robust regression in the sense that the result is not so affected by the presence 
of outliers. To do this, following Numerical Recipes (Press et al. 2007), we first perform a 
regression based on an assumption that the residuals follow a Lorentzian distribution rather 
than a Gaussian distribution. Then we derive a tentative relation and evaluate the intrinsic 
dispersion, which can be used to identify outliers. Eventually, removing the outliers, we perform 
ordinary chi square regression to derive a final relation with confidence intervals of fitting 
parameters. Here it should be noted that the robust regression does not give the confidence 
intervals, which is why we have to remove outliers and perform ordinary chi square regression. 

Our method can be summarized as the following five steps. 

1. Choose samples with small measurement uncertainties. 

2. Fit a relation using a robust regression based on a Lorentzian distribution of the scatter 
(§-3.1). 

3. Calculate the robust standard deviation of the residuals and estimate the intrinsic disper- 
sion of the relation (§-3.2). 

4. Identify outliers (§-3.3). 

5. Remove the outliers and perform ordinary chi square regression on the remaining samples. 

We find that the first step is important for the data analysis of GRB whose measurement 
uncertainties are not uniform. The largest uncertainty about E p among our gold samples is 
aE P /E p ttl, while the smallest one is aE p /E p « 0.01. Even if there would be different populations 
in the diagram of GRBs, large observational uncertainties would make it difficult to distinguish 
them. 

The goal of this section is to describe each step of our method to derive the E p -Tl-L p 
relation and estimate the intrinsic dispersion. We will show a demonstration of our method by 
performing Monte Carlo simulations in Appendix 1. 

3.1. Robust regression 

Before explaining robust regression, let us review ordinary least-square regression. We 
assume a linear correlation between E p , Tl and L p in logarithmic scale as, 

\ogL p (E p ,T L ) = A + B\og (E p /E p ) + Clog (T L /f L ) , (2) 

where A, B and C are the parameters of the model. We adjust these parameters to maximize 
the likelihood function given by, 
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N j 

P(A 1 B 1 C)<xY[ex V (--z*)AL, (3) 
log L Pji - log L p (E p>i ,T Lt i) ^ 



where AL is an arbitrary small constant and a- lxA is the intrinsic dispersion of the correlation. 
The 2C factor in front of <J\ og L vi comes from the fact that the definition of Tl includes L p . 
Maximizing this likelihood function is equivalent to minimizing chi square function, 

N 

X 2 (AB,C)=J2zt. (5) 
i=i 

Because not only the numerator but also denominator of z$ depend on model parameters, this 
chi square function is not a linear function of model parameters. Therefore nonlinear regression 
algorithm, Levenberg-Marquardt Method (Levenberg 1944; Marquart 1963), is used to find the 
best fit parameters. If a int is not known in advance or there are extra components of error which 
cause the scatter, <7 int is adjusted to hold x^ in /d.o.f. = 1. However we should point out that this 
procedure to estimate <7i nt is based on an assumption that the intrinsic dispersion around the 
model follows a Gaussian distribution. Therefore, contamination of outliers and/or different 
populations will make the estimation incorrect. We thus develop a more reasonable way to 
estimate the intrinsic dispersion of the correlation, which does not depend on the existence of 
outliers and/or different populations (§-3.2). 

Let us move to robust regression. In this paper, we assume the scatter of the samples 
around the relation follows a Lorentzian distribution, rather than a Gaussian as suggested by 
Numerical Recipes (Press et al. 2007). In this case, the likelihood function is, 

jv 1 

P(A,B,C)<xH-— T ^AL. (6) 

i=l l '2 Z i 

This distribution has wide tails so that the existence of outliers is common. We maximize this 
likelihood function, or equivalently, minimize the sum of the negative logarithms. Ignoring 
constant terms, this means minimization of the Lorentzian merit function defined by, 



M{A,B,C)=YM 



N r x 



= 1 



1 + -z? 



(7) 



Using this merit function with a int = 0, a tentative set of model parameters are determined. 
Regression based on a Lorentzian distribution is much more robust than one based on a 
Gaussian distribution, because a Lorentzian distribution has wide tails and the contribution of 
outliers to the merit function is highly suppressed. 

Although the robust regression can find reasonable best-fit values of the model param- 
eters, it doesn't provide reliable confidence intervals for them, that is, it can't be used to 
compare the fit of different sets of parameters. This is why we perform an ordinary least-square 
regression in the final step. 
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3.2. Estimation of intrinsic dispersion 

To identify outliers, we need to measure how much each sample deviates from the re- 
lation, which is evaluated by Zi in Eq. (4). The value of the intrinsic dispersion a- int is then 
necessary and we estimate it as follows. 

If the residuals follow a Gaussian distribution, the distribution of zi follows a Gaussian 
distribution whose mean is zero and standard deviation is unity. In this case, we can expect 
that the standard deviation is estimated as, 

N N 

£{log£ P) i - togL p (J5 Pii) T Li4 )} 2 =X;{(l + 2Cy? ogLpi + B 2 al gEpi + CPo^J+Ntfa 

i=l i=l 

However, in the current case, the l.h.s. of Eq. (8) is strongly influenced by outliers. Thus, 
we replace it with a robust standard deviation of the residuals (ctrsd), for which we adopt the 
median absolute deviation (MAD), 

_ median[|logL p>i -logL p (£'p ii ,T Lii )|] 

aRSD = 0^745 • (9) 

Here the factor 0.6745 comes from the fact that 50% of a Gaussian distribution lies in 0.6745 
standard deviation of the mean. Obviously this orsd is not affected so much by the presence 
of outliers. It should be noted that, if the scatter around the relation follows a Gaussian 
distribution, this ctrsd is equivalent to the normal standard deviation. 

If we identify the l.h.s. of Eq. (8) with Acr RSD , we can estimate the af nt as follows, 

A' 



2 2 
a mt = a RSD 



1 

- £ {(1 + 2C)af osLp i + B^ Epi + C 2 *l gT J. (10) 

We use this value of <r int to calculate and then detect outliers. In Appendix 1, we demonstrate 
that this way of calculating the intrinsic dispersion is reasonable. 

3.3. Detecting outliers 

Let us explain how to detect outliers. Because we have obtained a tentative set of model 
parameters and the intrinsic dispersion, we can compute, for each sample, t = \zi\ and the 
two-tailed P-value from the t distribution with (N — K) degrees of freedom, where N and K 
are the numbers of samples and model parameters, respectively. Following Motulsky & Brown 
2006, we adopt the False Discovery Rate (FDR) method of (Benjamini & Hochberg 1995). We 
need to determine the threshold value of P-value which decides whether a sample is an outlier 
or not. To do this, Motulsky & Brown 2006 ranked the P-value from high to low and defined 
the threshold value of P-value for i-th sample as a, = Q(N — i + l)/N where Q is an arbitrary 
number less than unity. This means that the i-th sample with Pj < ctj is regarded as an outlier. 
In this paper, we simply define the threshold value as, 

a = Q/N, (11) 

for all i, which makes our criterion more conservative, and we take Q = 0.1. For this choice of 
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Q, we mistakenly regard 10 /N% of samples as outliers on average even if there are actually no 
outliers. 

Here, we have estimated the intrinsic dispersion for data set including outliers and then 
it is slightly larger than the value without outliers. For more reasonable estimation of the 
intrinsic dispersion and outlier detection, we repeat the steps 2-4 again. 

By now, we have shown a new method to derive correlation in presence of many outliers 
and/or different populations. We apply this method to estimate the best fit parameters and 
the intrinsic dispersion of the E p -Tl,-L p relation. 

4. The improved E p —T L —L p plane 

In this section we apply the method which we have developed in the previous section to 
the data described in section 2. Here, we take a time resolution of L p as a free parameter of 
the relation, and find the most favored time resolution of the relation. 

We use only the platinum data set to derive the E p -T L -L p relation as discussed in 
section 3. In the left of Fig. 1, we show the minimum of the Lorentzian merit function as a 
function of time resolution of L p . We found the most favored time resolution is approximately 
3 seconds with cx int = 0, and six outliers (080319B, 081222, 090328, 090926, 091003, 091127) are 
identified. After removing these outliers, we perform chi square regression with <7i nt = 0. Then 
we obtained following result, 

/ P \ l-90±0.036 / T \ -0.12±0.053 

J _ -I n 52.50±0.017 ( ffr \ ( J-L \ /..„>. 

p ~ vl0 2 - 656 keVy VlO°- 95 seJ ' 1 } 

with x^ in /dof = 12.16/9. The most favored time resolution is 2.51q;2- m Fig. 2, we show the 
best-fit E p -Tl-L p relation of (18-6) platinum GRBs. The red points indicate the GRBs which 
are used to estimate parameters, and the green points are the GRBs which are eliminated as 
outliers. We should emphasize that, because our method removes outliers in automatic way, 
there is no artificial choice of outliers. 

If we include the four mid-outliers (080319B, 081222, 090926, 091127), the chi square 
value becomes unity when <7; nt = 0.20. As discussed in section 3 and appendix 1, if we assume 
the intrinsic dispersion cr int = 0.20, the chance probability to obtain the value cr int = 0.0 using 
our method with 15 samples is approximately 1% (See top central of figure 6). Thus, we can 
conclude that these four GRBs surely outliers of the E p -Tl-L p plane, or at least the true 
intrinsic dispersion of the E p -T^-L p plane is much less than 0.20. 

In top of figure 3, we showed not only the data used for calibration but also other 
data (gold and bronze data sets). Bottom left figure indicates the histogram of weighting 
residuals (zj) and bottom right indicates the histogram of unweighting residuals (log L p i — 
log L p (E P j,TL,i)). As one can see, the gold data are consistent with best-fit model within 
statistical error, but the most of bronze data distribute below the best-fit line of Eq. (12). It 
might be explained by the fact that the peak energies estimated by the CPL model becomes 
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Fig. 1. The minimum Lorcntzian merit function (left) and chi square function (right) as a function of 
time resolution of L p . It should be noted that the Lorentzian merit function is calculated using all 18 
platinum samples, while the chi square function is calculated after eliminating 6 outliers (080319B, 081222, 
090328, 090926, 091003, 091127). 

systematically larger than that estimated by the Band model (Kaneko et al. 2006; Krimm et al. 
2009). 

5. Summary and Discussion 

In this paper, we considered two possible origins of the systematic error in the Ep-T^- 
L p relation. The first is concerned with the calculation of the spectral peak energies (E p ) and 
we removed samples whose E p was estimated by CPL model. The second is concerned with 
the definition of the peak luminosity (L p ) which were conventionally estimated by 1-second 
in observer time. We converted the time resolution to the time scale in GRB rest frame. 
Furthermore, because there might be some outliers like 980425 and some short GRBs, we 
developed a new method using robust regression and following outliers rejection. Using our 
method, we found that the E p -T L -L p relation is the tightest around 3-second peak luminosity 
with <7 int = 0. 

We briefly discuss the properties of six outliers detected by our method. First, we should 
note that three of them (090328, 090926, 091003) are detected by Fermi /LAT. Especially GRB 
090328 has long-lived GeV emission after the prompt emission and GRB 090926 has extra high 
energy component. It might be possible that Fermi /LAT selectively observe some kind of 
GRBs which are harder than normal GRBs. Second, GRB 0803 19B is one of the remarkable 
GRBs ever observed because of the optical emission bright enough to be seen by naked-eyes. 
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Fig. 2. The E p -Tl~L p diagram for 18 platinum samples. The red points indicate GRBs which is used to 
derive the relation and the green points indicate 6 outliers eliminated by our method described in section 3. 
Solid line indicates the best-fit model in Eq. (12). 

Although the others (081222, 091127) seem to be normal GRBs, these might belong to peculiar 
GRBs mentioned above. 

Of course, because the number of the sample is small, our result must be confirmed 
by future experiments. The collaboration of Swift and Fermi can provide a few GRBs with 
measured redshift and small E p error with <je p /E p < 0.1. Future mission SVOM will be able 
to determine the redshift and E p by itself, and can confirm our result with a large number 
of data. In future studies, we should study the classification of GRBs using comprehensive 
observation of both prompt and afterglow emission. The Ep-T^-Lp plane might become useful 
discriminator of GRBs like HR diagram of stars. Once we establish the relation, GRBs might 
be a very powerful and unique distance indicator to probe the high redshift universe. 
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Fig. 3. (Top): Same as figure 2, but gold data (blue points) and bronze data (light blue points) are also 
plotted, (left bottom): The histogram of the weighting residual from the best-fit line. The color follows 
the top figure, (right bottom): The histogram of the unweighting residual from the best-fit line. 
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Table 1. Intrinsic property of GRBs for 18 platinum data set'' 
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b The data with E p fitted with the Band function and aE p /E p < 0.1 
f 2.496 second peak luminosity in GRB rest frame 
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Table 2. Intrinsic property of GRBs for 14 gold data set 
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b The data with E p fitted with the Band function and ue p /E p > 0.1 
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Table 3. Intrinsic property of GRBs for 54 bronze data set'' 
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Table 3. (Continued.) 
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b The data with E p fitted with the CPL function 
f 2.496 second peak luminosity in GRB rest frame 
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Appendix 1. Monte Carlo simulation 



Here we show the validity and limitation of our method using Monte Carlo simulations. 
For simplicity, we consider a relation between 2 quantities, X and Y, in this section. 

A. 1.1. A case with two populations with different normalizations 

Here we consider two populations whose relations are, respectively, 



logF 1 = 48 + 21ogX 1 , (Al) 

logF 2 = 48.3 + 21ogX 2 (A2) 

with o"i n t = for both of them. 

We generate mock data according to the following equations, 

= C/(100, 1500) (A3) 

a Xll /X h = a Yii /Y h = U (<7 min , a max ) (1 £ % x £ N x ) (A4) 

Y h =48 + 2 logX 4l + G(0, v/^iogy^+^iogx^+^t) (A5) 

X i2 = Z7(100, 1500) (A6) 

a Xl2 /X i2 = a Yl2 /Y i2 = U (a min , a max ) (1 £ i 2 £ N 2 ) (A7) 

Y l2 = 48.3 + 21ogX 2 + G(0, ^+5^+0 (A8) 



where C/ (min, max) represents a random number from a uniform distribution between min 
and max, and G(m,o"sD) represents a random number from a Gaussian distribution which has 
mean value m and standard deviation (Jsd- We generate N x and X 2 (N x > N 2 ) samples for each 
population with 100(7^% and 100cr^% observational errors in X and Y, respectively. Then, 
following the method described in §-3, we obtain a tentative set of best-fit model parameters 
and the intrinsic dispersion, and identify outliers. We assume the fitting function of the form 
logF = A + B\og(X/X). 

First, we consider with relatively small observational errors, (min, max) = 

(0.01,0.05). In figure 4, we show the result of 1000 simulations for (N U N 2 ) = (12,3) (top), 
(40,10) (middle) and (80,20) (bottom), respectively. The left figures show the result of one re- 
alization and each line represents the best-fit line for ordinary regression (dashed line), robust 
regression (dotted line), and our method (dash-dotted line). The assumed relation Eq (Al) 
is also indicated by the thick line. All lines except the ordinary regression are overlapping. 
The white points indicate the outliers which are detected by our method. The central figures 
show the histogram of the estimated <7i nt . As the number of sample increases, the value of a- mt 
converges to the true value. We show distribution of the value of the estimated parameters 
in the right figures. The parameters estimated after outlier elimination are indicated by black 
points, while the ones obtained by the normal regression with all samples are indicated by gray 
points. Apparently our method is more effective than ordinary regression. We summarized the 
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Table 4. Results of Monte Carlo simulations with two populations. 



-^total 




Cint (fiducial) 


iV out (fiducial) 


15 


[0.01,0.05] 


0.02 (0) 


2.95 (3) 


50 


[0.01,0.05] 


0.01 (0) 


10.09 (10) 


100 


[0.01,0.05] 


0.01 (0) 


20.1 (20) 



Table 5. Results of Monte Carlo simulations with two populations with large observational uncertainties or large fraction 
of outliers. 



AWi [CTmm,c max ] a ixA (fiducial) N out (fiducial) 
100 [0.1,0.2] 0.09 (0) 0.62 (20) 

100 [0.01,0.05] 0.01 (0) 26.8 (40) 



parameters and results in table 4. Figure 4 and table 4 show that our method gives more rea- 
sonable results than ordinary regression in a case with two populations whose normalizations 
are slightly different like Cepheid variables. 

Next, we consider a case with relatively large observational errors, (min,max) = (0.1,0.2), 
to show the limitation of our method. Top left of figure 5 shows that there are no points 
eliminated as outlier, although the assumed number of outliers is 20. Thus, if observational 
uncertainties are larger than the difference in the normalization of the two relations, it is difficult 
to detect outliers. This is why we use only samples with small observational uncertainties. 
Likewise, if there are too many outliers (Bottom), 40% rather than 20%, it is also difficult to 
detect outliers correctly. The result is summarized in table 5. 

A. 1.2. A case with intrinsic dispersion 

Next we consider a single population with the correlation defined as, 

logy = 48 + 2 logX. (A9) 

Here we assume that this relation has the intrinsic dispersion of cr int = 0.20, which is much 
larger than observational errors (< 0.05). 

We generate mock data according to the following equations, 

X i = [7(100, 1500) (A10) 
o- Xi /Xi = cryjYi = £/(cT min ,CT max ) (All) 

Y t = 48 + 21ogX, + G(0, V / a 1 2 ogyi + J BV 1 2 ogXi + cr mt ). (A12) 

In figure 6, we show the result for iV= 15 (top), 50 (middle), 100 (bottom), respectively 
The meaning of points and lines are the same as in figure 4. For these simulations, there 
are little points eliminated as outliers. As the number of sample increases, the value of cr int 
converges to the fiducial value. It is indicated that our method return reasonable values even 
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Fig. 4. Results of Monte Carlo simulations with two populations whose normalization factors are slightly 
different. The number of samples is (iVi,iV 2 ) = (12,3) (top), (40,10) (middle) and (80,20) (bottom), 
respectively. The left figures show the result of one realization and each line represents the best-fit line 
for ordinary regression (dashed line), robust regression (dotted line), and our method (dash-dotted line). 
The assumed relation Eq (Al) is also indicated by the thick line. All lines except the ordinary regression 
are overlapping. The white points indicate the outliers which are detected by our method. The central 
figures show the histogram of the estimated a m t- As the number of sample increases, the value of dint 
converges to the true value. We show distribution of the value of the estimated parameters in the right 
figures. The parameters estimated after outlier elimination are indicated by black points, while the ones 
obtained by the normal regression with all samples are indicated by gray points. See also table 4. 
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Fig. 5. Results of Monte Carlo simulations with two populations whose normalization factors are slightly 
different. The assumed two correlations are the same as in figure 4, but we consider two cases where 
observational uncertainties are larger (top) and where the fraction of outliers is larger (bottom). The 
meaning of points and lines are the same as in figure 4. See also table 5. 

Table 6. Monte Carlo simulation with intrinsic dispersion. 



iVtotal 




(Tint (fiducial) 


iV ou t (fiducial) 


15 


[0.01,0.05] 


0.20 (0.20) 


1.18 (0) 


50 


[0.01,0.05] 


0.20 (0.20) 


0.51 (0) 


100 


[0.01,0.05] 


0.20 (0.20) 


0.28 (0) 



if the correlation has intrinsic dispersion which follows Gaussian distribution. We summarized 
the parameters and result in table 6. 

A. 1.3. A case with two populations with different intrinsic dispersions 
Next we consider two populations with the same correlation, 

log ^,=48 + 2 \ogXj (j = l,2), (A13) 

20 





Fig. 6. The results of Monte Carlo simulation with the intrinsic dispersion <j m t = 0.2 for N = 15 (top), 
50 (middle), 100 (bottom), respectively. The meaning of points and lines are the same as in figure 4. For 
these simulations, there are little points eliminated as outliers. As the number of samples increases, the 
value of dint converge to the fiducial value. See also table 6. 
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Table 7. Results of Monte Carlo simulations with uniform distribution outliers. 



-Wtotal 


[^"min i *2"max] 


cri n t (fiducial) 


N out (fiducial) 


15 


[0.01,0.05] 


0.032 (0) 


2.7 (2) 


50 


[0.01,0.05] 


0.011 (0) 


9.15 (10) 


100 


[0.01,0.05] 


0.011 (0) 


18.13 (20) 



but with different intrinsic dispersion. For j = 1, we assume a- mt = and, for j = 2, we assume 



the scatter around the relation is uniform between —1 and 1. 

We generate mock data according to the following equations, 

= [7(100, 1500) (A14) 

ax n /X h = a Yil /Y h = U (a mm , <r max ) (1 < % x < N ± ) (A15) 

Y h = 48 + 21ogX n + G(0, v/<y n + ^ 2 <x n + °L) (A16) 

X i2 = [7(100,1500) (A17) 

a Xi2 /X i2 = a Yi2 /Y i2 = U (a min , a max ) (1 = i 2 < N 2 ) (A18) 

F i2 =48 + 21ogX i2 + f/(-l,l) (A19) 



In figure 7, we show the result for N — 15 (top), 50 (middle), 100 (bottom), respectively. 
The meaning of points and lines are the same as in figure 4. We summarized the parameters 
and result in table 7. 

These result indicate that our method return reasonable values even if there are two 
types of correlation whose intrinsic dispersions are different. 
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